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ABSTRACT 

We study the stability of poloidal magnetic fields anchored in a thin accretion disc. The 
two-dimensional hydrodynamics in the disc plane is followed by a grid-based numeri- 
cal simulation including the vertically integrated magnetic forces. The 3-dimensional 
magnetic field outside the disc is calculated in a potential field approximation from the 
magnetic flux density distribution in the disc. For uniformly rotating discs we confirm 
numerically the existence of the interchange instability as predicted by Spruit, Stehle 
& Papaloizou (1995). In agreement with predictions from the shearing sheet model, 
discs with Keplerian rotation are found to be stabilized by the shear, as long as the 
contribution of magnetic forces to support against gravity is small. When this support 
becomes significant, we find a global instability which transports angular momentum 
outward and allows mass to accrete inward. The instability takes the form of a to = 1 
rotating 'crescent', reminiscent of the purely hydrodynamic nonlinear instability pre- 
viously found in pressure-supported discs. A model where the initial surface mass 
density S(r) and B^ir) decrease with radius as power laws shows transient mass ac- 
cretion during about 6 orbital periods, and settles into a state with surface density and 
field strength decreasing approximately exponentially with radius. We argue that this 
instability is likely to be the main angular momentum transport mechanism in discs 
with a poloidal magnetic field sufhciently strong to suppress magnetic turbulence. It 
may be especially relevant in jet-producing discs. 
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1 INTRODUCTION 

In the formation and evolution of Young Stellar Objects 
magnetic fields are thought to play a key role. Magnetic 
fields, beside rotation and thermal pressure, stabilize molec- 
ular clouds from gravitational infall. By magnetic braking, 
thermal cooling or diffusive magnetic processes like ambipo- 
lar diffusion a supercritical cloud can form. The subsequent 
collapse, which proceeds preferentially along the rotation 
axis and the magnetic field lines, results in a dense central 
object, i.e. a protostar and an accretion disc (Mestel 1965; 
Shu et al. 1993; see Tomisaka 1995 for a numerical study). 
After the first, dynamical collapse we expect the disc around 
the central object to be still threaded by a fraction of the 
interstellar magnetic flux of the original cloud. 

Though microscopic diffusion processes are important 
for the evolution of magnetic fields in star forming regions, 
effective diffusion of matter across field lines is also possible 
through instabilities driven by the magnetic field itself. In 



the molecular cloud cores where the magnetic field is thought 
to disengage itself from the matter, the gravitational and 
magnetic energy densities are believed to be of comparable 
magnitude (see McKee et al. 1993 for a review). This is just 
the condition under which magnetic instabilities can operate 
at rates competitive with the gravitational collapse rate. It 
would be somewhat of a coincidence if microscopic diffusion 
were always the dominant process, and effective diffusion 
by magnetically-driven instabilities never played a role in 
the entire contraction process from cloud to star. This is 
one important reason to study the possible effects of non- 
axisymmetric magnetic instabilities in the process of star 
formation. 

Another reason is provided by the possible connection 
of magnetic fields in discs with outflows and jets. A remnant 
of the cloud core's magnetic flux, anchored in the accretion 
disc, represents a poloidal magnetic field component with 
the right geometry to launch and accelerate disc gas along 
the magnetic field lines (Bisnovatyi-Kogan and Ruzmaikin 
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1976; Blandford & Payne 1982). P Cygni line profiles in the 
spectra of Young Stellar Objects might find their explana- 
tion in such disc winds (Edwards, Ray & Mundt 1993). 

In addition to their role in accelerating outflows, 
poloidal magnetic fields anchored in the accretion disc may 
also play an important role in the collimation of bipolar out- 
fiows and jets (Blandford 1993; Spruit, Fogglizio & Stehle 
1997). We refer to Hughes (1991) for an extensive introduc- 
tion into the observations and physics of beams and jets. 
For a review and tutorial introduction to the magnetic ac- 
celeration model see Spruit 1996; for some issues of current 
interest see Ogilvie and Livio (1998), Cao and Spruit (1994, 
2000). 

The ability of the disc to produce a magnetically accel- 
erated outfiow depends rather critically on the strength and 
radial distribution of the poloidal field at its surface. This 
distribution, in turn, is determined by the rate at which an 
advected large scale magnetic field in the disc is able to dif- 
fuse outward against the inward drift velocity compressing 
it. 

Van Ballegooijen (1989) studied the radial transport of 
magnetic field lines assuming an isotropic turbulent viscos- 
ity u, related to an equally isotropic magnetic diffusivity r) 
by a constant ratio, the magnetic Prandtl number. Along 
the same lines Lubow, Papaloizou & Pringle (1994) studied 
the magnetic field dragging by turbulent processes in an ac- 
cretion disc which was initially threaded by an externally 
imposed magnetic field. The conclusion from these studies 
is that inward dragging of an external field is difficult to 
achieve, if the assumption of an isotropic magnetic Prandtl 
number of order unity holds. It follows that one should ex- 
pect to see, in this case, only internally generated fields like 
those obtained in numerical simulations (Stone et al. 1996; 
Brandenburg et al. 1995). 

If the initial magnetic field is sufficiently strong to con- 
tribute to support against gravity, however, as appears to be 
the case in cloud cores, it is likely to suppress these dynamo 
processes, since a weak-field instability (Balbus & Hawley 
1992) is thought to be an essential ingredient in this kind of 
dynamo. In this case, which we shall call here the strong field 
case, the most plausible source of turbulence in the disc are 
non-axisymmetric instabilities caused by the strong mag- 
netic field itself. 

One well-known instability is interchange, a local insta- 
bility driven by the magnetic field energy, and which plays an 
important role in controlled fusion devices and in solar mag- 
netic fields. The example of solar prominences shows that 
these instabilities can be quite efficient in transporting mass 
across field lines (cf. Priest 1982). Linear interchange insta- 
bility has been studied previously by Spruit & Taam (1990) 
for uniformly rotating discs, and for differentially rotating 
discs by Spruit, Stehle & Papaloizou (1996) and Lubow & 
Spruit (1996), in a shearing sheet approximation. The insta- 
bility appears when the ratio Bz/S of the vertical magnetic 
field strength to the surface mass density decreases with dis- 
tance from the centre of the disc. Its behavior is similar to 
convection, with the gradient of Bz/S replacing the entropy 
gradient. 

As in the case of solar prominences, it is the magnetic 
tension force due to the curvature of the field lines in the r — 
2-plane that drives the instability in geometrically thin discs 
(Anzer 1967; Spruit and Taam 1992; for a detailed analysis 



of the magnetic forces in thin discs see Ogilvie 1997). Since 
non-axisymmetry is also crucial, this makes the magnetic 
field three-dimensional. Our study thus differs conceptually 
from the numerical study by Kaisig Tajima and Lovelace 
(1992), where a 2-dimensional cylindrical configuration with 
vertical magnetic field lines was assumed. In that case, the 
only magnetic force is the magnetic pressure gradient. 

A full three-dimensional numerical treatment of the 
problem is made difficult by the very large range of char- 
acteristic speeds expected in the problem. Inside the disc, 
the Alvfen speed is not larger than the sound speed, but 
in the low density regions outside the disc it can easily ap- 
proach the speed of light. We show here how this difficulty 
can be turned into an advantage, such that only a part of 
the problem needs to be done in 3 dimensions, while the 
numerically most demanding parts can be done using only 
the two dimensions in the disc plane. 

In the above, we have introduced the magnetic disc 
model with the example of protostellar discs. The physics 
addressed by our calculations, however, is equally applica- 
ble to discs in compact binaries or active galactic nuclei. 



1.1 poloidal fields vs. dynamo-generated fields 

The magnetic problem we study is in several ways comple- 
mentary to that of the simulations done by Stone et al. 
(1996) or Brandenburg et al. (1995), and Hawley (2000). 
There the magnetic fields are generated locally by a dynamo 
process, which can already start from a very weak intitial 
field. If a strong field is present initially, such that the mag- 
netic energy density is comparable to the thermal energy 
density, the dynamo process is suppressed. This limit can 
be written in the form 



47rS 



(1) 



where E is the surface mass density and Cs the sound speed. 

A small scale dynamo process like the magnetic turbu- 
lence seen in these simulations does not create a net fiux 
of field lines through the disc (Hawley 2000), and the over- 
all field structure is therefore at least of quadrupole order 
at large distance from the disc. For advected magnetic field 
lines, as they are advected radially and anchored in the disc, 
the total magnetic flux through the disc is non-zero. The 
global structure of the advected magnetic field, as seen from 
large distances, is thus close to a dipole magnetic field dis- 
tribution. The radial force exerted by such a field is pre- 
dominantly the tension force, whose magnitude integrated 
over the disc thickness is BrBz/2-K. A poloidal field can in 
principle exist up to strengths such that this radial force 
starts contributing signiflcantly to the support against grav- 
ity. This limit can be written as 



47rS 



^ fi r-. 



(2) 



This limit on B^ is a factor flr/cs larger than the strength 
(^) at which dynamo-generated fields are suppressed. Once 
strong poloidal fields exist in a disc, they suppress the mag- 
netic turbulence that would make them diffuse out of the 
disc by van Ballegooijen's argument. There is thus a large 
range in parameter space where a poloidal field in a disc 
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would be affected only by its own internal instabilities. These 
are the subject of the present investigation. 



2 TWO-DIMENSIONAL DISCS WITH 
THREE-DIMENSIONAL 
MAGNETOSPHERES 

We neglect the self-gravity of the disc so that the central 
star is the only cause of gravitational force acting on the 
disc. Following Spruit & Taam (1990, hereafter ST, see also 
Tagger et al. 1990 for a similar discussion) we assume a geo- 
metrically thin disc with an internal velocity field that does 
not depend on the vertical coordinate (perpendicular to the 
disc). The equations of motion can then be integrated across 
the disc, resulting in a two-dimensional problem confined to 
the plane of the disc. This plane can in principle be of ar- 
bitrary time-dependent shape, as in ST, but we limit the 
calculations here to the case of a plane disc, without dis- 
placements of the disc surface in the vertical direction (i.e. 
without corrugations or 'bending modes'. For results on such 
modes, in the same approximation for the magnetic field, see 
Agapitou et al. 1997). 

One may wonder when the approximation of height- 
independent velocity fields in the disc is justified. The 
assumption clearly eliminates the possibility of dynamo- 
generation of magnetic fields, but is appropriate for strong 
fields in the sense discussed in the preceding section, namely 
fields which contribute, more than the pressure force, to sup- 
port against gravity. For such strong fields, winding-up of 
field lines inside the disc can be ignored, because the small 
amount of differential rotation encountered by the field line 
as it crosses from one side of the disc to the other is not 
enough to bend field lines significantly. Instead, the differen- 
tial rotation across the disc will adjust, due to the magnetic 
forces, such that the rotation rate is constant along a field 
line. In the thin-disc limit, the amount of differential rotation 
across the disc vanishes as H/r, justifying the assumption 
made (for a discussion of the thin limit of magnetic discs, 
from a more mathematical point of view, see Ogilvie 1997) . 

Because of the thin disc limit taken, the dominant mag- 
netic force in the problem is the curvature force due to the 
bend of the field lines crossing the disc (see ST). Outside 
the disc, the magnetic forces dominate over fiuid forces on 
account of the low gas density. We call this region of low 
plasma-/3 the 'magnetosphere' of the disc (not to be confused 
with the magnetosphere of an accreting magnetic star). We 
simplify the physics by assuming the field in this magneto- 
sphere to be current free so that it is derivable from a mag- 
netic potential ip. This is justified if the matter density in 
the magnetosphere is sufficiently low. The magnetospheric 
field in then approxmately force free. The case of the solar 
corona shows that such a field is in practice also close to a 
potential (current free) field. This is due to the fact that at 
low P, 'forced' processes (e.g. Parker 1979, for recent numer- 
ical simulations see Galsgaard and Nordlund 1997) are fsist, 
and keep the degree of twisting in the field low. 

In this approximation, all currents are confined to the 
plane of the disc, and are proportional to the jump in the 
tangential field components across the disc. These tangen- 
tial components {B^, Br) are determined by the normal field 
component through the solution of the (3-dimensional) po- 



tential problem in the magnetosphere, for which the normal 
component B^ provides the boundary condition. 

The magnetic forces acting on the fluid in the disc plane 
are given by the difference in the magnetic stress acting on 
the upper and lower surfaces of the disc. They are propor- 
tional to the product of the tangential and normal compo- 
nents of the magnetic field at the disc surface (see ST for 
details). 

In the computations, closed inner and outer boundaries 
are used for the disc (vr — 0, so that no mass or magnetic 
flux crosses these boundaries. Thus the total magnetic flux 
through the disc and the total disc mass are constant in 
time. 

Apart from the addition of the magnetic force term, the 
hydrodynamical problem is the same as in ordinary two- 
dimensional disc hydrodynamics. We use an Eulerian grid 
with the van Leer (1977) scheme for upwind differencing. 
An outline of the model and its basic assumptions has been 
given previously in Stehle (1997). 



3 EQUATIONS 

To describe the non-axisymmetric disc response resulting 
from large scale magnetic fields, we adopt a thin-disc ap- 
proximation in which the vertical velocity vanishes. We use 
a cylindrical coordinate system (r, (f), z) defined with the ori- 
gin at the central mass M. The mass surface-density is 
E = pdz with p = p(r, 4>, z) the gas density. 

In the thin disc limit H/r <ti 1 where H is the disc 
thickness the only contribution of the Lorentz force per unit 
surface area is due the magnetic tension of the field lines 
(ST). This is because for magnetic fields varying on a length 
scale L 3> -ff, the curvature force is larger than the magnetic 
pressure gradient by a factor L/H. Then the radial compo- 
nent of the equation of motion for the disc, assumed to be 
inviscid, reads 

"57 + ^■'"5""' "5T (3) 

ot or r o4> r 

1 dP B4B,] 

T, dr 47rS ^' 

where and = rfl are the radial and azimuthal compo- 
nents of the velocity vector and gris the acceleration of grav- 
ity due to the central star. B^ is the vertical magnetic field 
component at disc midplane and [Br] = B^ — B~ = 2B^ 
the jump of the field vector from above (B^) to below (-Br) 
the disc plane. Neglecting any possible disc warps it is as- 
sumed that the field geometry is antisymmetric with respect 
to the disc plane (i.e. B^ = B^ and B^ — —B^, see ST). 
For use in what follows we define the magnetic acceleration 
in radial direction, 

= B,B+/27rE. (4) 

The azimuthal equation of motion is 

dv4, dv^ dvtj, ViV^ 

"oT + ^'■"a ' ^ 

ot or r 0(p r 

1 dP B,[B^] 

^ Urdcj) 47rE 

The continuity equation is 
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dt r dr 



d(t)' 



(6) 



A similar equation holds for the vertical component of the 
magnetic field in the limit of complete flux freezing: 



dB. 1 d 

ot r or 



1 d 



(B^rv^) = 0, 



(7) 



which expresses that the magnetic flux density Bz is con- 
served. 

The gas pressure is computed from the vertically inte- 
grated internal energy e assuming an ideal gas for the equa- 
tion of state, P = (7 — l)e, with ratio of specific heats 7. 
For the computations reported here, the value 7 = 1.4 was 
used. The adiabatic evolution of the internal energy is given 
by 



de 19, , 1 a , , 

at r or r-' oq) 



(8) 



To close our set of equations we have to determine the 
inclination of the magnetic field lines to the surface of the 
accretion disc, i.e. we need to derive B^{r,(f>) and B'^{r,(f)). 
Their values determine the magnetic forces in the equations 
of motion. By our assumption of a potential field in the mag- 
netosphere, electric currents only exist within the plane of 
the accretion disc, and the vertical component of the electric 
current vanishes everywhere. 

The current is therefore of the form 



jk{r,4>,z) =jk{r,(l>)6{z) and j^{r,4>,z) = 



(9) 



where S{z) is the Dirac-delta function. A subscript h denotes 
vectors parallel to the plane of the accretion disc. Disk winds 
and ionized particles in the disc magnetosphere, neglected in 
our model, will certainly contribute to current fiows in the 
disc magnetosphere. The calculation of the 3-dimensional 
disc magnetosphere in force-free approximation is, however, 
beyond the present computational feasibility. The assump- 
tion of a potential field for the structure of the disc magne- 
tosphere is equivalent to assuming that the magnetosphere 
is sufficiently close to its minimum energy state. 

In the thin disc limit, the magnetic field at the disc 
surface {B^ , Sj, Bz) is connected to the disc currents by 



R+-2^,- R+- 2^,- 

- — J0, B^ - > 



and 



(10) 



" n M!:Mz^dr'd<^' (11) 

^ri„ ^ 



where the second term on the right hand side (RHS) of eq. 
( prj ) is the current flow at the disc-boundaries and R is given 
by 



R =r + r' — 2rr' cos{<j) — <j)'). 



(12) 



Eq. ( |ll| ) is derived from a partial integration of the vec- 
tor potential fleld A in Coulomb gauge so that B = V x A, 
and (Landau & Lifshitz 1975): 



A{r, 



Ft 



(13) 



Given the magnetic flux distribution in the disc Bz{r, (f>) 
eq. (jl^) has to be inverted to yield the currents in the disc. 
The inversion is unique by applying the fact that currents 
do not accumulate, divj — 0: 

-!-('■>) + -^j0=O. (14) 
r Or r Oq> 

This connects the radial jr and azimuthal compo- 
nent of the current, and has to be solved together with the 
inversion of eq. ([ll[). 

As our accretion disc has a central hole, the solution of 
the potential field problem with Bz{r-,(j)) given in the disc 
is no longer unique, since the domain space is multiply con- 
nected. To any solution of the inversion problem (fnl), an 
arbitrary multiple of the solution for Bz = G can be added. 
This special solution consists of closed field lines wrapping 
through the hole of the disc and around its outer edge with- 
out crossing the disc itself (like the windings of a ring-core 
coil). We thus have additionally to specify the number of 
field lines which pass through the central hole, i.e. the total 
magnetic fiux through the hole of the disc (see also Lubow, 
Papaloizou & Pringle 1994) 



r Bz{r, <l),z — 0) dr d(/!>. 



(15) 



A Fourier Transform of eq. (|15|) shows that only the axisym- 
metric component contributes to the magnetic flux through 
the disc hole whereas all other components cancel. The de- 
gree of freedom introduced by the presence of a hole thus 
enters only into the computation of the axisymmetric com- 
ponent of the field. 

Eqs. (^, (^) and ^ specify the time evolution of E, B^ 
and e when the velocities Vr and in the plane of the accre- 
tion disc are known. These follow from Euler's equation (^ 
and (^) which includes thermal, gravitational and magnetic 
forces. To solve for the magnetic forces we invert equation 
( pr] ) with the differential constraint (|l^) for the unknown 
currents, and an assumed value for the magnetic fiux though 
the hole. Eq. (|l^) gives the relation between B^ , and j,/,, 
jr and thus the magnetic forces are determined. In the next 
chapter we show how we solve the hydrodynamic equations 
and the magnetospheric potential problem numerically. 



4 NUMERICAL METHOD 

We solve the equations on an Eulerian grid with equidistant 
spacing in the r and coordinates. The inner disc rim is 
at rin = 0.1 Tout with rout the radius of the outer edge of 
the disc. The number of grid points and grid spacing in 
radial direction are Ur and Ar = (rout — rin)/nr. We use 
a staggered grid such that the scalar quantities S, Sz, e are 
defined at the cell centres and the vector quantities (u,-, v^), 
{Bi,B^) and {jti,,ji) at the cell boundaries. The equations 
are used in dimensionless form, as follows. As unit of length 
we take the outer disc radius, as the unit of time the inverse 
of the Keplerian rotation rate at the outer edge. Thus, the 
Keplerian velocity at the outer edge is unity, and the time 
for one Keplerian orbital period of the outer edge is T = 2t7. 
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For all the calculations presented here we use rir x = 
156 X 128 grid cells in r and direction respectively. We also 
performed some calculations on a smaller grid of 64 x 64 grid 
cells. Comparing these models to calculations done with a 
higher resolution we find only differences on the order of the 
applied grid spacing. We are thus convinced that the models 
presented here are resolved sufficiently. 



the inner and outer edges of the disc, i.e. Vr = dP/dr = 0. 
These are sufficient for the present calculations in which 
only the short-term evolution of the disc is followed. For 
the longer-term evolution, one would want to use conditions 
that allow accretion to take place through the boundaries. 
This is beyond the scope of the present study. 



4.1 The hydrodynamic part 

The hydrodynamic part of the calculations is done with a 
natural extension to the scheme described in Stehle & Spruit 
(1999). 

The equations are written in conservative form. Terms 
in the equations are divided into advection and source terms. 
The advection from one grid cell to the other is done with 
the upwind-differencing scheme of van Leer (1977). Thermal 
pressure, gravitational forces and compressional heating are 
calculated following Stone & Norman (1992). Viscous pro- 
cesses and radiative cooling from the surface of the accretion 
disc are neglected. 

The induction equation, i.e. eq. (^, is an additional 
equation compared with the hydrodynamic case. It has the 
same forms as the continuity equation, and is treated nu- 
merically in the same way. 

The magnetic force is a new term in the equations of 
motion. Contrary to the pressure force, which is calculated 
by a local derivative, the magnetic forces are given by the 
bending of the magnetic fields where they pass through the 
accretion disc. The inclination of the magnetic field lines to 
the accretion disc depends, however, on the global magnetic 
field structure in the magnetosphere and only when this is 
known, the magnetic forces be calculated. We introduce the 
numerical method to calculate the magnetic forces in Section 

(0)- 

The time step is controlled by the Courant-Friedrich- 
Levy condition. In addition to the sound speed, there is a 
magnetic wave speed in the problem. The wave is compres- 
sive, but since the restoring force results from the change 
in the external potential upon compression of the field lines 
rather than of the magnetic pressure itself, it is a dispersive 
wave. The phase velocity is (ST, Tagger et al. 1990) 

- (16) 



k 



\/27rEfc' 



where k — 2n/X is the wavenumber and A the wavelength. 
The group speed, which carries the wave information, is 
a factor 2 lower. The highest magnetic wave speed in the 
discretized equations is therefore obtained for the highest 
wavenumber that can be represented by the grid. By the 
Nyquist theorem, this is fc = n/Ar. The magnetic wave 
speed Vm that enters the Courant condition is thus 

= ^ . (17) 

27r^2E/Ar 

A Courant factor of 0.75 was found to be sufficient for nu- 
merical stability. In most cases, however, we find the time 
step to be controlled by the azimuthal velocity v^{rin) at 
the inner accretion disc rim. The additional magnetic wave 
constrains the time step only for models where the magnetic 
acceleration Qm approaches the acceleration of gravity g. 
The boundary conditions used are solid boundaries at 



4.2 The solver for the disc magnetosphere 

The magnetic forces are calculated at each time step from 
the magnetic fiux distribution Bz{r,(f)) in the disc. This in- 
volves two steps. First, eq. (^ij) is inverted to obtain the cur- 
rents jh from the flux distribution B^. In the second step, the 
forces BzB^ /2tv, B^B'^ /2tv are computed using eqs. (p^, 
and added to the hydrodynamic forces. 

The inversion of eq. ( pT[ ) has two complications. Be- 
cause of charge conservation the current is not an arbitrary 
function of r and (f), but must satisfy eq. (p^), i.e. divj — 0. 
This condition can be used to eliminate one of the current 
components jV and in favor of the other. Eq. (^l|) can 
then be read as an integral equation determining one of the 
current components in terms of the vertical field component. 

Second, the region on which the computations are done 
is not simply connected. The central hole in the computa- 
tional domain generates an additional parameter in the po- 
tential problem, namely the net magnetic fiux through the 
hole ^f. For the present exploratory calculations, we use fixed 
boundary conditions. Hence Vr = at the boundaries, and 
with the induction equation (^ no magnetic fiux enters or 
leaves through the edges of the disc. This also implies that 
is constant in time. More general boundary conditions 
are possible that take account of the advection of magnetic 
fiux into the hole. These will be needed in more realistic 
calculations of discs accreting on to protostars, for example. 



4-2.1 Fourier decomposition in (f> 

Since the equations for the potential problem are homoge- 
neous in the azimuth 4>, Fourier transforms can be used in 
this direction. Since they are also linear, the Fourier com- 
ponents do not mix, and one can solve for each of the 
Fourier components separately. If the number of azimuthal 
grid points or Fourier modes is n^, this reduces the com- 
puting effort required by a factor of the order compared 
with straightforward discretization in (j). By using Fourier 
decomposition, the computing effort for the potential prob- 
lem can be kept at a level comparable to the hydrodynamic 
parts of the calculation. 

Thus we write the magnetic flux distribution as 

B.(r, 0) =B°^Y^ (Sr-' sin(m0) + B^'" cos(m0)) , (18) 



where B°{r), B^'''{r) and Bl^'^ir) are only functions of r. 
Similar equations hold for the currents jr, j^. 
The Fourier amplitudes are given by 



B™'''(r) = - / B^{r-,(j))cos{m' 
Jo 



(19) 



© 2001 RAS, MNRAS 000, 



6 R. Stehle, H.C. Spruit 



and similarly for B^'^. Using this becomes 

dr'[r'Mr',(P')]-d^,jrir',cl>') 



(20) 



cos{m4 



(^2 + r'2 - 2rr' cos((/) - (p'))^^^ 
This can be written as 



where 



1 

TTCr 



dr'd<^' 



(21) 







cosijrul 



(22) 



(1 + a;2 - 2j:cos((?i'))^''^' 

We evaluate this function numerically. For m = it can be 
expressed in terms of the complete elliptical function of the 
second kind F{x) (Gradstein & Ryzhik 1981), 



Ko{x) = 4:xF{x) with x = niin(a:, l/s). 



(23) 



Substituting the Fourier expansions of jr and and 
integrating over 0', we get 



[9r'(r'j7' 



tjr'1^m(r'A)dr'. (24) 



For m 7^ the charge conservation condition ( |l^ ) can be 
used to eliminate j^, which yields 



(njr 



1 

cr 



dr' 



(25) 



-a.,(r'9,,r'jr-=) - m9^,jr'= K^^r' /r) 



where T'"''^(r,ji.) is an abbreviation for the integral operator. 
A similar equation holds for the _B™'° component. 
The axially symmetric component reads 

("■out 



T\r,U) = -^j d,,{r'jl)Ko{r'/T)AT'. (26) 
The flux through the hole in the middle of the disc is 

= 2:*(j^) = 27r [ rB°dr = (27) 
Jo 

^ ~c / 3r'(r'j^)A'o(r7r-)dr'dr. 



Eqs. (|2^) and ( |2q ) are integral equations for the (n^ — 1) 
Fourier amplitudes j7^''^{r) and jl"'^{r) of the radial current, 
and the azimuthal current distribution j'^{r). To solve these, 
we take flnite differences in r, which turns each of the eqs. 
( p5[ ) and (j26|) into a set of linear algebraic equations. The 
matrices involved are flxed in time and need to be inverted 
only once. If the number of grid points in r is rir, the comput- 
ing effort for the potential problem is therefore of the order 
n^n^ per time step, or ~ rir per grid point and time step. 
Since the number of operations per grid point and time step 
in the hydrodynamic part of the calculation is a substantial, 
but fixed number independent of rii, the computing expense 
for the potential problem does not dominate the overall ex- 
pense except for very large numbers of radial grid points. In 
fact, for a grid of rir xn^ = 156 x 128 grid cells the magnetic 
part of the calculation takes about 50% of the CPU-time. 



4-2.2 Discretization in r 

The discretization in r of the integrals in (|2^) fc ( pij ) is 
different for the non-axisymmetric components (Ga) and the 
axisymmetric component (p^. The Fourier amplitudes _B™ 
of the field are naturally defined at the same radial positions 
as the values of the field B^ itself, i.e. at the centres of 
the cells defined in the hydrodynamic part of the scheme. To 
discretize (|6|) , the currents are defined at the boundaries 
''i+i/2 = {fi + ''i+i)/2 between these cells. Because of the 
topology of the domain used (a disc with a hole) , there is one 
more of such boundaries than there are cell centres. Since 
the boundary conditions do not impose constraints on the 
azimuthal current at the boundaries, there is then also one 
more current than there are field amplitudes B^ . This 
additional degree of freedom is balanced by the the hole- 
flux condition (|l5|), which arises from the same topological 
property. 



4-2.3 The axisymmetric component 

To evaluate the integral in (En) we interpolate the current 
linearly between the values ..^at its grid points: 

j^.^(r) = (x-l)jO + xJ^'+i, (28) 

where x = {r - rj+1/2) / {rj+i/2 -rj+i/2)- Inserting this into 
(pij), B1 is a linear function of the currents jj", with coeffi- 
cients B%: 



(29) 



Evaluation of these coefficients involves integrals of the type 



I (r) = / {a + Pr')Ko{r'/ri)dr'. 



(30) 



The elliptic function involved in Ko is evaluated by poly- 
nomial approximation (Abramowitz and Stegun 1984), the 
integral by Bulirsch's algorithm (Press et al. 1995). 

Similarly, assigning the index i = to the hole-flux "if, 
( p?! ) can be written as 



(31) 



The coefficients B^^ then form a square matrix of dimen- 
sion Ur + 1, relating a vector of magnetic variables, b" — 
(*,B°(ri),...B^(r„^)), to the currents: 



bi — ^ ^ Bij Jj . 



(32) 



Since we use a fixed Eulerian grid, the matrix elements 
Bfj are fixed and the inversion of the matrix can be done 
once and for all for a given computational grid. This inver- 
sion is done by the LU decomposition method (Press et al. 
1995). 



4.2.4 The non-axisymmetric components 

The procedure for the non-axisymmetric components is very 
similar to that for the axisymmetric component, except that 
does not appear because it is already determined by the 
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5 UNIFORMLY ROTATING DISCS 
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Figure 1. An example showing the accuracy of the magneto- 
spheric field calculation for an m = 1 mode. The assumed az- 
imuthal field at the top surface of the disc B^ \{r) (dashed) agrees 
with the reconstructed values (solid) to 1 % for this grid with 100 
points in r. 



axisymmetric component. For m 7^ the currents J™, like 
the , are defined at the cell centres, so that there is an 
equal number of each. 

Because we have used a partial integration in deriving 
(p5|), and have used charge conservation to eliminate j^, the 
integrand in (|2^) contains a second derivative with respect 
to r. In order to evaluate it at the same order of accuracy as 
the axisymmetric coefiicients, a third order interpolation is 
needed. We choose cubic spline interpolation between neigh- 
bouring grid-cells. The coefficients then involve expressions 
of the type 



jr' — 



'it'm(x)dx 



(33) 



with I — 1,2,3. Coefiicients up to m = 64 were evaluated 
with an accuracy of 10~^. 

4-2.5 Test calculations 

To test the accuracy of our numerical solution of the po- 
tential field problem, we apply the inversion of eq. (^l|) to 
a known field configuration. We proceed by first specifying 
Ji(r) analytically. Then we derive Sz(r) from a numerical 
integration of eq. ([ll|). The integration yields B^(rs) at the 
grid cell centres ri with an accuracy of ~ 10"''. These val- 
ues are used to invert with our magnetosphere-solver to 
reconstruct the disc currents and field strengths 

which are then compared with the original function values 
at these points. 

An example is shown in Fig. (^) . This is an m = 1 mode 
with the current distribution 



3rA 



(r- 



■(^ou 



(34) 



with Tin = O.Olrout. The reconstructed current distribution 
matches the input distribution within an accuracy of 1%, 
for the 100-point grid used. 

Additional tests of accuracy of the code as a whole were 
done by comparing results at different resolutions. 



The linear stability of uniformly rotating discs with mag- 
netic fields of the type considered here was studied by ST. 
The condition for instability of the interchange type is 



_ S+B. d , ^ 



< 0. 



(35) 



Disks where a > everywhere are stable to interchange 
instability. The derivation of this condition does not take 
account of possible global instabilities. 

In this section we study two examples, one of a disc 
which is stable (model 1) and one that is unstable (model 
2) according to condition (jj^). 

Uniformly rotating discs are set up numerically by mod- 
ifying the gravitational potential <I>(r) such that magnetic, 
centrifugal and gravitational forces are just balanced for the 
case v^{r) =const. While this case is of limited astrophysical 
interest, it serves to test the agreement with the predictions 
from the linear theory, to check for possible global modes 
not covered by condition (^), and to get an impression of 
the nonlinear development of interchange instability in the 
present case. 

We study the discs in a frame of reference which coro- 
tates with the disc. The corresponding non-inertial terms 
are added to the equations in section 2 (see also Stehle & 
Spruit 1999). Both models presented here are advanced in 
time with zero magnetic flux through the central hole, i.e. 
with 'I' = const. = 0. 

The uniformly rotating case as defined above is a one- 
parameter family, governed by the ratio Cu of magnetic to 
the centrifugal forces: 



5m 

a 



f rBA 



(36) 



Other parameters such as the amplitudes of the central po- 
tential, the magnetic field strength and the surface density 
can be scaled out of the equations. In cases where Cu ^ 1 
magnetic forces are unimportant and the disc rotates freely. 
In the Cfl.SG Cu — ^ 1 magnetic forces dominate and rotation 
can be neglected. 

In the results shown, the gas pressure included in the 
calculations for numerical reasons (section ^) has only a 
small influence. 



5.1 Model 1: d(Bz/E) /dr = 

According to eq. ( psj ) discs with d(Bz/E) /dr = are ex- 
pected to be stable against the interchange instability. The 
initial density and magnetic field distributions are specified 
as Gaussian humps in r [E(r) ~ exp(— (r — 0.5 rout)^/A^)], 
with a maximum at r = 0.5 rout and width A = 0.1 rout. 

Models with four different field strength are followed: 
a weakly magnetized disc (cu — 910"*), and one with a 
high magnetic support (model la, Cu = 225). Intermediate 
cases are chosen with Cu = 0.09 and Cu = 9 . We perturb the 
initial stationary, axisymmetric models with a low amplitude 
(~ 10" '^w^ (rout)) point-to-point noise in v-^. 

We find the extreme cases Cu ^ and Cu ^ cxa to 
be stable. The total kinetic energy in the radial velocity 
component, (i.e. £kin,r = / Sw^dF) was constant for the 
whole calculation of ~15 disc orbits. 
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Figure 2. The time evolution of tlie Icinetic energy -Bjjin r for 
tile uniformly rotating model sequence. The duration of one disc 
revolution at rout is 27rT. The interchange instability grows on a 
dynamical time scale determined by the magnetic field strength. 
Both the growth rate and the saturation amplitude increase with 
the field strength. 



In the intermediate cases Cu ~ 1 a very weak form in- 
stability was observed, with characteristics different from 
an interchange. The energy in the radial motions i?kin.r in- 
creased during the first 15 orbits by a factor ~ 100-1000 . 
A global disc pattern is excited, showing spiral arms with 
m ~ 20. The wave saturates at a strength (B^jB^ ~ 10~®- 
IQ-*. The waves are of low amplitude and do not measur- 
ably transport angular momentum. After a linear rise time 
of some orbits the disc is found to be stationary again, 
i.e. 9tS(r, 0) = 0, even though it is now slightly non- 
axisymmetric. 

The exact nature origin of this weak instability is not 
entirely clear at the moment. In any case, the amplitude 
of the motions observed is low compared to those of the 
instabilities described below, and are not relevant for actual 
accretion discs where the magnetic support is always less 
than gravity, Cu < 1. 

5.2 Model 2: d (B^/E) /dr < 

Next we study accretion discs where condition ( ^sj ) predicts 
the interchange instability to be present. 

We choose S =const. and B^ir) decreasing with r as a 
power law 

B4r) = B,,o{ ■ (37) 

V rout / 

With this choice, the instability parameter a has a minimum 
at r « 0.2 Tout- We expect the instability to first manifest 
itself near this radius. 

As in model 1 we perturb the initial axisymmetric model 
with a low amplitude point-to-point noise in Vi. 

Instability is found in all models of this sequence. Fig. 
(^) shows the time evolution of i5kin,r = / ^Sw^ dF for these 
models. The model parameters for this sequence are shown 
in Table H. At first -Ekin.r increases exponentially. -Ekin.r sat- 



Table 1. Parameters for model sequence 2: initial values of the 
degree of instability am = min(a) and the ratio of magnetic to 
centrifugal acceleration Cu at r = 0.2 Tout, and growth rate 7e of 
the kinetic energy -Bkin r- 



nr. 


dm 


Cu 


7e 


2a 


-1.10 


1.1 


16.4 


2b 


-0.36 


0.5 


9.6 


2c 


-0.16 


0.2 


6.5 


2d 


-0.05 


0.05 


3.4 



urates after several orbits and the subsequent evolution is 
highly non-linear. 

The instability causes a significant redistribution of the 
disc mass. This is illustrated in Fig. (^, which shows the 
evolution of the surface mass density for model 2c. It is seen 
that the instability first operates at r ~ 0.2 rout as expected 
from the local minimum of a. A high mode number m ~ 15 
dominates at first. The pattern of motions resembles that 
of convective cells or the plumes of a Rayleigh- Taylor insta- 
bility. The influence of rotation is evident in the asymmetry 
of the plumes. With time the instability is seen to spread 
over the whole disc. The small instability cells merge and 
grow in size as they drift to larger radii, as is characteristic 
of Rayleigh- Taylor instability. 

We identify this instability with the interchange insta- 
bility as the disc pattern looks very similar to convective 
cells as predicted by Spruit et al. (1995), and starts at the 
point where the linear analysis predicts the disc to be most 
unstable. 

We conclude that the instability in the uniformly rotat- 
ing case sets in as predicted from linear theory and has the 
nonlinear development of an interchange instability. 



6 NON-UNIFORMLY ROTATING DISCS 

We now study accretion discs revolving around the gravita- 
tional field of a point mass. According to the linear analysis 
of Spruit et al. (1995) and Lubow & Spruit (1995), shear 
due to differential rotation acts as a stabilizing factor on the 
interchange instability. This analysis predicts that instabil- 
ity appears only in regions of the disc where magnetic forces 
contribute significantly to support against gravity. The pre- 
dicted linear growth is algebraic (a power law of time) rather 
than exponential. 

We study the evolution of two different initial setups. 
First we follow a model where the magnetic field decreases 
as a power of radius (model 3) and then a case where it de- 
creases exponentially (model 4). The initial structure of the 
models is summarized in Tab. (^) and The initial v^{r) 
is found from the radial force balance between magnetic, 
gravitational and centrifugal forces. We then perturb Vi by 
point-to-point low amplitude noise and follow the subse- 
quent disc evolution numerically. As before all models are 
calculated with zero magnetic flux through the central hole 
of the disc, * = 0. 
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Figure 3. Interchange instability in a uniformly rotating disc with -Bz/S decreasing outward, and an initially uniform surface density 
S (model 2b). The instability starts near to the inner edge of the accretion disc, where the magnetic instability parameter a is largest, 
and spreads over the whole disc in a few orbits. 
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Figure 4. The initial axisymmetric condition of model 3a. Initial 
surface mass (a) magnetic field (b) and azimuthal velocity distri- 
bution (c). The ratio of the magnetic ((/m) and gravitational (g) 
accelerations is shown in (d). 



c 




Figure 5. Growth of the instability for model sequence, showing 
evolution of the power Pm in the first 4 Fourier modes, relative to 
the power in the axisymmetric mode. Curves are shifted vertically 
by factors of 10 since the saturation levels are nearly the same. 
The power in the 4 modes is comparable, except during the hump 
at t/T ~ 40 in model 3a ('outburst') when m = 1 dominates. 



6.1 Model 3: 



-5/4 and E 



-3/2 



For the initial state in model sequence 3we choose a sur- 
face density varying as E ~ r~^^^. In order to contain the 
effects of the instability within the computational domain 
as much as possible, we choose a magnetic field distribution 
with strength vanishing towards the boundaries (Fig. ^). 

The degree of support against gravity by the magnetic 
field, as measured by the ratio pm/ggj, increases with radius 
and has a maximum at 0.85 rout (Fig. Md). The parameters of 



Table 2. Parameters for model sequence 3. The relative support 
of the disc by magnetic forces (3m/Sg)max, the initial growth rate 
7Ekin,ri the ratio TZ of the total magnetic flux through the disc to 
the total mass, and the mass accretion rates in units of the total 
disc mass Af^jsc per time unit T. 



model nr. 


(9m/5)max 


7Ekin,r 




M(r < 0.3rout) 


3a (t/T < 35) 


0.32 


27.1 


0.93 


2.0 lO-i Mdisc/T 


3a (t/T ~ 40) 








1.710-3 Mdisc/T 


3b 


0.14 


6.9 


0.63 


1.210-5Mdi,c/r 


3c 


0.035 


0.96 


0.31 


~410-«M'disc/r 



the model sequence are summarized in Tab. ()^. The param- 
eter TZ specifies the mean magnetic flux per surface-mass: 



7^: 



/dis 



BzdF 



/dis 



Edi^ 



(38) 



and can be used to compare model calculations. 

Instead of interchange type instability, we find, in all 
three cases that the initial set up is unstable to a global, non- 
axisymmetric instability. The wave pattern of the instability 
can be traced from one edge of the disc to the other (see Fig. 
1^. Initially the kinetic energy i?kiii,r growth exponentially 
on a dynamical time-scale. The initial growth rate, as given 
by 7Ekin,r = d lu iJkin.r /dt , is largest for the disc with the 
highest magnetic support. 

Fig. (@) shows the time evolution of the power 



/;°"' B2 rdr 



(39) 



in Fourier mode m of the field strength B^, integrated over 
the whole disc, and relative to the axisymmetric component 
m = 0. In all three calculations we find that the modes m = 
1-4 grow equally fast. The relative power saturates nearly 
at the same level, independent of the degree of magnetic 
support. The relative power in the first 4 modes is similar. 
For the calculation with the highest magnetic field strength, 
however, (model 3a) the m = 1 component dominates for a 
period of about 15 orbits around t/T ~ 40 (the numerical 
time urut T = l/Slout, where flout is the orbital frequency 
at the outer edge). 

Fig. ^ shows the evolution of the mass in the inner 
disc (rin < r < 0.3 rout) in units of the total disc mass 
Afdisc- Mass piles up in the inner part of the disc, the faster 
the higher the magnetic support. This is accompanied by an 
outward transport of angular momentum by the magnetic 
instability. For model 3a we find a roughly linear increase of 
the inner disc mass with time corresponding to an accretion 
rate Mdisc(r < Q.3rout) ~ 2 10~*MdiBc/T, but superposed on 
this trend is a much more active 'outburst' around t ~ 40. 
During this active episode the accretion rate is about 10 
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Figure 6. The evolution of the mass in the inner disc (rin < 
r < 0.3 Tout) for models of sequence 3. Only model 3a shows 
a significant mass accretion towards the central star. It shows 
a period of enhanced accretion around t = 40 during which a 
prominent m = 1 spiral arm is present, cf. Figures 5 and 7 



times higher (i.e. for 35 < t/T < 40). Fig. (|) shows that 
during the outburst Pi is larger by at least a factor 3-5 
compared to the other modes, and by a factor ~ 10 larger 
than in the preceding phase. Figs. show that a strong 
m=l spiral wave, traveling outward from a crescent-shaped 
disturbance in the inner disc, is present during the outburst. 

The time-scale for mass accretion in model 3a is longer, 
by a factor lO'^-lO*, than the dynamical time-scale. For discs 
with less magnetic support the mass accretion time-scales 
are so long that we have been unable to follow their evolution 
beyond the initial development of the instabilities. 

The evolution of the disc pattern as seen in Bz is shown 
in Fig. (^) . The corresponding images of the surface density 
E are found in Stehle (1997). It is seen that the global insta- 
bility starts with rather high mode numbers, m ~ 5. ..8 (Fig. 
1^). Subsequently the waves are wound up (Fig. |^) and it 
is only later that the modes m = 1...4 become dominant. 
After about 30 orbits of the outer disc edge, a prominent 
m — 1 spiral arm develops. It causes mass accretion rates 
10 times higher than during the preceding phase. The rel- 
ative strength in the m —1 component is largest near the 
inner disc edge. At r ~ 0.2 rout a prominent crescent-shaped 
field strength enhancement shows up (best seen in Fig. Mc). 
It is accompanied by a similar enhancement in the surface 
density (Stehle 1997). The crescent rotates approximately 
with the local orbital rate. An m = 1 wave travels outward 
from this rotating crescent. The maximum of the crescent 
moves in to smaller radii (compare Figs. ^ and e), where 
the rotation rates are higher. The outward traveling wave 
becomes correspondingly more tightly wound. 

This behavior is very reminiscent of a purely hydrody- 
namic form of global instability observed in hot, partially 
pressure-supported discs (Blaes and Hawley 1988; Rozyczka 
& Spruit 1993). This instability takes place only when the 
degree of support against gravity by pressure becomes no- 
ticeable, and it also takes the form of a crescent rotating 
at the local orbital rate. It generates shock waves travel- 
ing outwards and inwards. These waves cause the mass in 



the disc to spread, while the angular momentum lost by the 
crescent causes it to spiral in towards the centre, behaving 
much like a solid object in doing so. Though the waves in 
the present calculation are rather different from hydrody- 
namic shock waves, we suspect that the same mechanism is 
at work. The peculiar behavior of the crescent mode and the 
fact that it appears only at certain phases suggests that it 
is a basically nonlinear phenomenon, not related directly to 
the linear global modes of the system. The nature of this 
phenomenon warrants further study. 

With time, mass piles up near to the inner and outer 
edge of the accretion disc. The effect of the instability is thus 
much like that of viscous spreading, but it must be stressed 
that this is due to a global transport of angular momentum 
by the spiral wave, which can not be reduced to the action of 
a local viscosity. After the approximately 10 outer disc orbits 
during which the m = 1 component dominated, its relative 
strength compared to the other modes decreases again and 
the phase of high mass accretion rates is finished (Fig. ^) . 

At the end of the outburst the disc density distribution 
has completely changed. It now decreases approximately ex- 
ponentially with radius rather than as a power law. The 
same is true for the distribution (Fig. |^ . 

The dissipation taking place during the redistribution 
of mass and magnetic flux causes the disc to heat up, so 
that the gas pressure increases over its initial low value. At 
the end of the calculation the Mach number of the orbital 
motion in the inner disc has decreased to values of 5-10. 

For the calculations done with magnetic support less 
than in model 3a the mass accretion rates are too low for 
significant redistribution of mass to occur over the 50 orbits 
we were able to follow. It is thus unclear if discs with less 
magnetic support also show outbursts, or if a threshold in 
the magnetic field strength exists below which the outburst 
mechanism cannot operate. 

Since the outburst was a transient but obviously very 
effective transporter of mass and angular momentum, one 
wonders what caused its decline. After the outburst the de- 
gree of support of the disc against gravity by the magnetic 
forces has decreased substantially. Since the amplitude of all 
global waves seen in our results increases sharply with the 
degree of magnetic support, it is possible that the outburst 
declined because the crescent instability operates only at 
sufficiently high degrees of support. This view is consistent 
with the properties of the purely hydrodynamic crescent in- 
stability. 

Another possibility is that an exponential dependence 
on radius is perhaps a more stable configuration, towards 
which the instability tends to develop. To test this possi- 
bility, we investigate in the next model a sequence of discs 
where E and decrease exponentially with radius. 

6.2 Model 4: B^ ~ exp(-r) and E ~ exp(-r) 

The final disc structure of model 3a motivates us to study 
the stability of accretion discs where the initial magnetic 
field and surface mass density decrease exponentially with 
radius. In this sequence we study two such cases, which differ 
only in the initial field strength. In model 4a we choose a 
field with magnetic support ((?m/(?)max = 0.36 at maximum 
and in model 4b one with (3m/ff)max = 0.16. The initial 
structure of model 4a is shown in Fig. (|9|) . 
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Figure 7. Snapshots of the evolution of model 3a. Instability sets in with rather high m. Later the modes m =1-4 are strongest until at 
t/T ~ 35 the m = 1 component becomes dominant. During the presence of the m = 1 mode the mass accretion toward the central object 
is enhanced. The outward traveling m = 1 spiral is generated by a rotating crescent-shaped enhancement which drifts inward with time. 



Figure 8. The distribution of Br/Bz and In S at the end of model calculation 3a {t/T = 62.4). The surface density decreases exponentially 
and the inclination of the magnetic field lines towards the outer disc edge decreases. 



Table 3. As Tab. 3, but for model sequence 4 



model nr. 


{9m/ 9g)max 


7Ekin,r 


n 


M{r < O.Srout) 


4a 


0.36 


5.8 


3.0 


8.710~®M/T 


4b 


0.16 


3.7 


2.0 


6.210-<'M/T 
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Figure 9. The initial surface density (a) magnetic field (b) and 
azimuthal velocity distribution (c) of model 4a. B2,{r) ~ S(r) ~ 
exp(— r). Panel (d) shows the initial ratio of magnetic to gravita- 
tional acceleration in the model. See text for further details. 

We take Bz/E = const, for the whole disc. E(r) and 
Bz{r) decline linearly to zero at the inner edge. The initial 
axisymmetric distribution is again perturbed in Vr by a low 
amplitude point-to-point noise. 

Values for the ratio TZ are listed in Tab. (P) for mod- 
els 3 and in Tab. for models 4. A comparison of these 
values shows that the total magnetic flux through the disc 
in models 4 is higher than in models 3, for the same de- 
gree of magnetic support. This difference is partly caused 
by the rapid decline of the magnetic field strength towards 
the disc edges in models 3 and partly from the fact that 
d(Bz/E)/dr > for the initial distribution in model 3. 

Fig. ( p^ ) shows the time evolution of the relative power 
Pm in the Fourier modes m =1-4, for models 4a and 4b. In 
both models the relative power Pm again increases initially 
on a short, dynamical, time-scale and the mode strength at 
which they saturate is similar to what we observe in model 
sequence 3. As expected from model 3a at times after the 
outburst, Pm=i is comparable to Pm=3 and Pm=4 but signif- 
icantly weaker than Pm=2- The m = 2 mode appears to be 
the dominant mode during most time of our integration. 




10 20 30 40 50 
t/T 

Figure 10. The relative power Pm{t/T) for the modes m =1—4 
in models 4. The values for model 4b are vertically shifted by 
multiplying Pm with 0.1. Pm=2 is most dominant and the Pni=i 
component is weaker during most time of the integration. 
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0.3 Tout) in models 4. The accretion time-scales are large com- 
pared to models 3. 



Even though TZ is significantly higher than in model 
3, the mass accretion rates at the inner part of the disc 
are small compared to model 3. The time-scale to clear the 
disc mass completely at these rates is now about 10"" outer 
edge orbits. This can be seen in Fig. (0), where we plot 
the evolution of the mass in the innermost part of the disc, 
M(r < 0.3 Tout) in units of the total disc mass Mdisc (see also 
Tab. ^). It is also interesting to note that the mass accretion 
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Figure 12. Radial velocity amplitudes at the end of model cal- 
culation 4a (t/T = 51.1). Unit of velocity is f^,Kepler(''out). The 
tightly wound spiral arms can be traced from one edge of the disc 
to the other. 
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Figure 13. The azimuthally averaged surface density < S(r) > 
at the end of the calculation t/T = 51.1 (full line), compared with 
the initial distribution (dashed line). Significant redistribution of 
mass has taken place only in the outer parts of the disc. 

is now much less dependent on the strength of the magnetic 
field. The low values for M(r < 0.3 rout) are accompanied 
by a low level of power in the m — 1 mode. This agrees 
with model 3, where the highest mass accretion rates are 
also found at a time where Pm=i is large. 

In Fig. ( p^ ) we show the radial velocities of model 4a at 
t/T — 51.1, the time where we stopped the integration, Fig. 
( [l3| ) shows the azimuthally averaged surface density < E > 
at that time. 

The pattern shows tightly wound spiral arms which can 
be traced from one edge of the accretion disc to the other. 
The radial velocities are rather small, of the order 10~^ of 
the orbital velocity at the outer edge of the disc. The mass 
redistribution has been strongest for r ^ 0.7 rout where it 
has taken only a few binary orbits. This is just the region 
where the magnetic support was strongest. In comparison, 
the mass redistribution in the inner regions is small. 

In both sequences 3 and 4 the redistribution of mass and 
magnetic flux appears to be closely related to the degree of 
support against gravity by the magnetic field configuration. 
It appears in regions where Qm/g exceeds 5-10 per cent. The 
comparison also shows that, unlike interchange instability, 
the global instability that causes this redistribution is not 
directly related to the flux-to-mass distribution B/E. 



7 CONCLUSIONS AND DISCUSSION 

We have studied, by numerical simulation, the stability 
of accretion discs threaded by strong large scale magnetic 
fields, assuming that the discs are geometrically thin. The 
simulations solve the MHD equations for a vertically aver- 
aged accretion disc in the (r, (/>)-plane. The disc magneto- 
sphere is calculated in potential field approximation, i.e. we 



treat the disc as a current sheet and assume the magnetic 
field outside the disc to be current free. This approximation 
allows us to follow the evolution of a 3-dimensional field con- 
figuration with the computational effort of a 2-dimensional 
simulation. 

Previous analytical studies by Spruit & Taam (1990), 
predicted that uniformly rotating discs would show an 
interchange-type instability, a local instability that appears 
at any location where the magnetic field decreases with 
radius more rapidly than the surface density E. To test 
this prediction, as well as the stability of the numerical 
method, we first computed discs with initially uniform rota- 
tion. The results agree with the analytic stability condition 
and growth rates. No additional, non-local, forms of insta- 
bility were found in the simulations of initially uniformly 
rotating discs. The nonlinear development of the instability 
agrees with that expected of an interchange instability like 
convection or Ray leigh- Taylor. 

Linear analysis (Spruit et al. 1995; Lubow & Spruit 
1995) predicts that differentially rotating discs with approx- 
imately Keplerian rotation are much more stable than uni- 
formly rotating discs. It predicts that instability of the in- 
terchange type occurs only when the local shear rate is less 
than the growth rate of the instability in a uniformly rotat- 
ing but otherwise identical disc. For smooth distributions 
of Bz and E with r this is equivalent to the condition that 
the magnetic forces contribute significantly to the support 
of the disc against gravity (Spruit et al. 1995). Such discs 
are strongly magnetized in the sense that va 3> flH, where 
va is the Alfen speed at the midplane of the disc and H 
the disc thickness. The importance of this prediction is that 
it suggests that even quite strong poloidal magnetic fields 
might still be stable in an accretion disc. 

To see if this prediction holds up in a full numerical 
simulation, we have done a sequence of calculations for discs 
with approximately Keplerian rotation, in which the mag- 
netic contribution to support against gravity ranges from 
a few % to 36%. At these values, the linear results pre- 
dict a weak form of interchange instability, with perturba- 
tions slowly growing as a power of t. In contrast with the 
uniformly rotating case, however, no evidence of this form 
of instability was found in the simulations. Instead, a new, 
global, exponentially growing form of instability appears. It 
takes the form of tightly wound spiral arms which can be 
traced from one edge of the accretion disc to the other. This 
instability was found both in cases where the magnetic field 
and the surface mass decrease with radius as a power law 
(model sequence 3) and where they decrease exponentially 
(model sequence 4). Its presence seems to depend primarily 
on the degree of magnetic support of the disc; its amplitude 
is a steep function of this quantity. 

The instability acts similar to viscous spreading in the 
sense that an outward angular momentum transport takes 
place which allows mass accretion from the outer to the inner 
regions of the disc. Mass accretion is of the order of 10~^- 
10"'* Afdisc/r where Mdisc is the total disc mass and T the 
Kepler time-scale 1/Q at the outer edge of the accretion disc. 

In the case with the highest degree of magnetic support 
an additional, more violent form of instability is observed. It 
takes place during a limited period of ~ 10 disc orbits. The 
mass accretion rate in this episode is 10 times higher than 
the long term average in the simulation. During this time the 
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disc is perturbed by a strong one-armed spiral wave excited 
by a density enhancement rotating at the local orbital fre- 
quency at the location of highest magnetic support. At the 
end of the outburst the surface density and the magnetic 
field distribution in the disc decline exponentially with disc 
radius and the subsequent mass accretion is less efficient. 

The time-scale for the redistribution can be rather large 
compared to the dynamical time-scale. Especially for discs 
with low magnetic support this time-scale is so large that, 
extrapolating the results from the simulation, the final con- 
figuration would be obtained only after some 1000-100000 
disc orbits. Only when the magnetic support of the disc is 
of the order of the gravitational or centrifugal forces, does 
the time-scale become small enough that the redistribution 
of the disc mass proceeds on time-scales of several 10 disc 
orbits. 

With the present results we can not establish whether 
the instabilities found can also lead to instability and redis- 
tribution of mass and angular momentum at lower degrees 
of magnetic support, or is limited to stronger fields. The dif- 
ference between these possibilities is of obvious importance 
for accretion on longer time-scales. If the instabilities found 
here generally require magnetic support exceeding a few per 
cent, quite strong poloidal fields might exist in the inner 
regions of accretion discs. 

In the protostellar context, such strong fields may also 
be relevant for the problem of binary formation (cf. Sigalotti 
and Klapp 2000) 
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